Self-assembly of artificial microtubules 
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Understanding the complex self-assembly of biomacromolecules is a major outstanding question. 
Microtubules are one example of a biopolymer that possesses characteristics quite distinct from 
standard synthetic polymers that are derived from its hierarchical structure. In order to understand 
how to design and build artificial polymers that possess features similar to those of microtubules, 
we have initially studied the self-assembly of model monomers into a tubule geometry. Our model 
monomer has a wedge shape with lateral and vertical binding sites that are designed to form tubules. 
We used molecular dynamics simulations to study the assembly process for a range of binding site 
interaction strengths. In addition to determining the optimal regime for obtaining tubules, we 
have calculated a diagram of the structures that form over a wide range of interaction strengths. 
Unexpectedly, we find that the helical tubules form, even though the monomer geometry is designed 
for nonhelical tubules. We present the detailed dynamics of the tubule self-assembly process and 
show that the interaction strengths must be in a limited range to allow rearrangement within clusters. 
We extended previous theoretical methods to treat our system and to calculate the boundaries 
between different structures in the diagram. 



I. INTRODUCTION 

Self-assembly of macromolecules or nanoparticles into 
polymers or superlattices is an important route to pro- 
duce higher-level structures with distinct properties for 
numerous applications.—^ Various assembled shapes can 
be formed, such as spheres, sheets and tubes, depend- 
ing on the structural features encoded in building blocks, 
their mutual interactions, and possible guiding from 
environment. In biological systems, self-assembly of 
biomolecules into extended structures is very common 
and provides crucial functions. One example is the self- 
assembly of microtubules (MTs) from a- and /3-tubulin 
heterodimers in cells. As key components of cytoskelc- 
ton, MTs play crucial roles in cell structure and move- 
ment, intracelluar protein transport, and cell division.-' 4 
Besides being an important biological macromolecule, 
MTs are special polymers that possess several properties 
distinct from standard synthetic polymers, e.g. tubu- 
lar structure, reconfigurable through depolymerization- 
polymerization cycles, and the substrate for motor pro- 
teins. MT assembly is a very complex process and thus 
difficult to characterize. Understanding the assembly of 
macromolecules that possess a subset of these properties 
would be remarkable and is our more practical goal. Our 
group is interested in developing polymers that are "arti- 
ficial" MTs, i.e. possess many characteristics of MTs, but 
are assembled from different monomers and constituent 
molecules^ In this paper we present results of an initial 
simulation study of the assembly of artificial MTs. We fo- 
cus on wedge-shaped monomers designed to self-assemble 
into tubular structures upon mutual binding and inves- 
tigate how the assembled structures depend on binding 
strengths. 

Tubulin monomers that make up MTs are about 4 nm 
in size£ They self-assemble into a helical tubular struc- 



ture when their concentration is above certain threshold 
or when there exist appropriate nucleation sites, such as 
7-tubulin ring complex in eucaryotic cells Each heli- 
cal turn typically contains 13 tubulin dimers in vivo, but 
the number can vary from 11 to 15 in vitro£ The diam- 
eter of MTs is around 25 nm, and their length can be as 
long as 25 /im. 

Another type of biological macromolecules that can 
self-assemble into tubes are surface layer (S-layer) pro- 
teins, which make up the cell wall of prokaryotic or- 
ganisms (bacteria and archaea). S-layer proteins typi- 
cally crystallize into two-dimensional arrays with vari- 
ous symmetries j™2 However, particular S-layer proteins 
under the right conditions can form open tubes with 
the diameter being controlled by adding or removing 
amino acids Various assemblies of S-layer proteins 
can serve as templates to grow other extended nanostruc- 
tures such as superlattices of metallic clusters . 11 ' 14 ' 15 

Tubular structures are also found to form in 
self-assembly of macromolecules that are usually 
amphiphilic^i^— The location of hydrophobic and hy- 
drophilic groups and other non-covalent interaction sites 
results in the formation of tubes. The key issue is 
to understand what features (shapes, chirality, binding 
strengths, etc.) of the molecular building blocks under 
certain formation conditions are required to yield tubes 
with desired structures (e.g., the tube diameter) and 
physical properties. However, there are so many vari- 
ables that it becomes difficult to explore the whole phase 
space experimentally. This is where molecular simula- 
tions might be most useful to identify the range of pa- 
rameters (particularly binding strengths) most appropri- 
ate to form tubes with controlled structures. 

An important, related example of self-assembly that 
for which some important aspects have been determined 
is the formation of a viral capsid from just one protein^ 
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Recent theoretical and simulation work have addressed 
the nature of this self-assembly and revealed some sim- 
plifying aspects of the assembly.— ~— The assembly of the 
protein monomer into any possible structure, not just 
into the capsid with icosahedral symmetry, has been ex- 
amined. While in the continuum limit, the free energy 
minimum favors icosahedral symmetry, for smaller cap- 
sids there are structural parameters that drive the spe- 
cific geometry. 28 Rapaport used MD simulations to study 
the formation of viral capsids from capsomers modeled as 
truncated pyramids with interaction sites on their lateral 
surfacesi 31 ' 33 Besides showing that these capsomers are 
able to self-assemble into capsids, his work found impor- 
tant general design aspects of the capsomer for assembly 
to occur. He found that there is a narrow range of inter- 
action strengths for the assembly of capsids. In addition, 
the simulations showed that reversibility along the assem- 
bly pathway is crucial for efficient production of complete 
shells by suppressing kinetic traps. 

Similar to Rapaport^ we use a 3D monomer with a 
wedge shape that promotes assembly into tubes. Repro- 
ducing all the properties of MTs is out of the question and 
given our interest in artificial systems, not our main goal. 
One initial purpose is to build a minimal model of the 
self-assembly of tubular structures. From the viral capsid 
results, we expect there to be a narrow range of interac- 
tions that will yield tubule self-assembly. Determining 
this range is essential before any further, more complex 
issues can be addressed. Our wedge monomers have lat- 
eral and vertical (also called "longitudinal" in MT liter- 
ature) binding sites on their surfaces, whose interaction 
strengths are independently varied. We have thus de- 
termined the optimal range of interaction strengths be- 
tween monomers that achieve tubules and studied other 
structures that form more generally at other interaction 
strengths. We observed a rich set of assembled struc- 
tures reminiscent of assembly of the coat protein into the 
tobacco mosaic virusi^i 

The assembly of wedge-shaped monomers is in the class 
of associating fluids and involves living polymerization 
under certain circumstances. In living polymerization, 
polymer chain can grow and shrink; ultimately there is 
an equilibrium between single monomers and a distri- 
bution of polymer lengths. The systems of associating 
fluids or those undergoing living polymerization can be 
modeled as an ensemble of particles allowing monomers 
to associate to form dimers, trimers, and long chains 
or clusters. The theoretical description of the assem- 
bly of these systems usually starts from either Wertheim 
theor y 35 i 36 which is essentially a perturbation theory, or 
Flory-Huggins theory^ Our model wedge, particularly 
when only two opposite sides are active, is in the class 
of polymer systems. In this paper we have developed 
a Flory-Huggins type theory that particularly treats the 
features of our model to calculate the structure diagram 
as a function of the interaction parameters. 

The remainder of this paper is organized as follows. 
In Sec. II, the simulation methods are described. Then 



in Sec. Ill a thermodynamic theory, along the line of 
Flory-Huggins lattice model of polymerization, is devel- 
oped for the self-assembly of anisotropic objects, which 
will be used later to explain the assembly behavior ob- 
served in MD simulations. The MD results are presented 
in Sec. IV. Finally, discussion and conclusions are in- 
cluded in Sec. V and VI. 



II. SIMULATION METHODS 

In order to have a monomer that will self-assemble into 
a tubule structure, the monomer is chosen to have a trun- 
cated wedge shape (Fig. [T]) . To make a wedge, we started 
with 27 particles in 3 x 3 x 3 simple cube. The distance 
between two neighboring sites is ler in the initial cube. 
The cube is deformed into a wedge shape such that 13 
of them will join to form a closed ring. The back layer 
is unchanged, while the front and middle layers are com- 
pressed such that the angle made by the two sides fits 13 
monomers in a ring. The particles (gray spheres in Fig. 
[1} interact purely repulsively through the Lennard- Jones 
(LJ) potential 

U(r) = 4e [(a/r) 12 - (a/rf - (a/r c ) 12 + {a/r c f} , (1) 

where r is the distance between the center of two parti- 
cles, and e is the unit of energy. The potential is trun- 
cated at r c = I.Oct to make the interaction purely repul- 
sive, which indicates that two monomers will repel each 
other when they get very close. 




(a) (b) (c) 



FIG. 1: (a) A wedge-shaped monomer; (b) An ideal ring 
formed by 13 monomers; (c) An ideal tube with 3 rings 
stacked. Each monomer consists of 35 sites, of which 27 par- 
ticles (gray) form its back-bone and the other 8 sites on the 
surfaces (in color) are the locations of the attractive interac- 
tion centers. Attractions only act between sites with the same 
color. 

The attractive component of the monomer- monomer 
interaction occurs at sites located on the wedge surface, 
specifically on pairs of sites on opposite faces (colored 
spheres in Fig. [1}. To control the orientation of two 
monomers that bind, two distinct attractive sites are nec- 
essary on a surface in order to break the symmetry of the 
surface-surface interaction. In Fig. Q] sites with the same 
color are attractive and sites with different colors have 
no interaction (effectively repulsive due to the other in- 
teractions). The red and blue sites on the lateral sides 
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promote the formation of rings. In order for the rings to 
stack to form a tubular structure, the two vertical (top 
and bottom) faces also have a pair of attractive inter- 
action sites (cyan, green). The attractive interaction is 
given by 



U(r) = -A 



1 



cos 



(2) 



where A is the strength of the potential and the range is 
given by r a . With this potential the strength and range 
of the attraction can be independently adjusted. In this 
work, we will only use r a = la. The value of A for 
the lateral and vertical surfaces will be independently 
varied. We will denote the potential strengths for the 
lateral surfaces to be Al and for the vertical surfaces to 
be A v . 

The dynamics of self-assembly was obtained by per- 
forming MD simulations using the LAMMPS simulation 
packagei 38 ! 39 Each wedge molecule is treated as a rigid 
body. The equations of motion were integrated using a 
velocity- Verlet algorithm with a time step 5t = 0.005r, 
where r = airn/t) 1 ^ 2 is the unit of time and m is the mass 
of one particle. The simulations to calculate the struc- 
ture diagram were run for 5 x 10 5 to 2.5 x 10 6 r. Unless 
otherwise noted, the simulations are for systems com- 
posed of 1000 wedges. To create an initial state without 
overlap of wedge monomers, the positions of the wedge 
monomers was taken from an equilibrated LJ system at 
density 0.2er -3 . Each LJ particle was replaced by a ran- 
domly oriented wedge particle and the system was scaled 
by the wedge particle size resulting in the volume frac- 
tion of wedges being 3.85%. The temperature of the 
system is kept at 1.0e/fcB, where fee is the Boltzmann 
constant, with a Langevin thermostat of damping rate 
l.Or -1 . Since the important energy scale is fee? 1 and 
since e = k^T, we will use k^,T as the energy unit. 



III. THERMODYNAMIC THEORY OF 
ARTIFICIAL MICROTUBULE ASSEMBLY 

To help understand the simulation results, we adapt a 
Flory-Huggins theory to describe the assembly behavior 
of wedge-shaped monomersj 28 ' 30 Suppose the simulation 
box is divided into a three dimensional lattice with M 
cells, each of which can only accommodate at most one 
monomer. The entropy of the system is clearly deter- 
mined by the number of ways to put N monomers into 
these M cells. As introduced in Sec. [TTl monomers simu- 
lated in this paper can form two types of bonds, lateral 
and vertical. The lateral bonding drives monomers to 
form closed rings, while the vertical bonding leads to the 
formation of filaments. Ideally, the appropriate combi- 
nation of these two types of bonds yields tubes. 

To get the critical bonding strength required for 
self-assembly, we first consider the simple case where 
monomers assemble into linear chains (e.g. Al = 0). 
Treating the linear case first allows us to follow much of 



the previous theoretical works i 30 ' 37 i 40 Denote the number 
of p-segment chains (consisting of p monomers) as n p and 
assume the maximum value of p is p max < N, the first 
constraint is from the conservation of total monomers 
and reads 



pn p . 



(3) 



The Helmholtz free energy of our system can be written 
as (see Appendix for a detailed derivation) 



M 



(4) 



where (p — l)Ag is the total energy gain to form a p- 
segment chain, which has p — 1 bonds, z is the coordi- 
nation number of the lattice, and j3 = X/k^T. This free 
energy has a form very similar to that derived by Zandi 
et al. for viral capsid assembly. 30 Minimizing the above 
free energy with the constraint Eq. leads to the dis- 
tribution of chains as 



n p = zM exp((3Ag) exp [—fip (Ag — fj,)] 



(5) 



where /i is a Lagrange multiplier and has a physical mean- 
ing of the chemical potential. Note that the law of mass 
action is implicit in the above expression of n p . Since 
N/M is the volume ratio of all monomers, fj, can be de- 
termined by combining the distribution n p and Eq. ([3]). 
Let x = exp [—/3 (Ag — jj)], we obtain 



z exp (fiAg) ^ P ;r 
P =i 



v — 



N 
M' 



(6) 



Complete the summation and take the limit p mayi — > oo, 
also notice that x < 1 to be physically meaningful, we 
arrive at the equation 



z exp (/3 Ag) 



(x iy 



N 



(7) 



from which x can be computed from known N/M and 
Ag. Then the chemical potential /i is given by /i = Ag + 
(3~ 1 lnx. With /j determined, the distribution of chains 
n p is easily computed. If we identify z _1 exp(— /3Ag) as 
the equilibrium constant of dimerization, Eq. (J7]) becomes 
identical to Eq. (5) in Ref. 40, which was derived directly 
from the law of mass action. 

The assembly is purely controlled by the bonding en- 
ergy Ag when the total number of starting monomers 
and the volume of the box are fixed. At small Ag, 
monomers dominate and ri\ is much larger than all n p 's 
with p > 2. With the increase of the magnitude of 
(< 0), n\ monotonically decreases. The onset of assem- 
bly is identified with the condition that half of the all 
starting monomers are in the assembled state, i.e., in p- 
segment chains with p > 2, while the remaining half is 
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in the form of free monomers*^ Using Eq. ([5]) and the 
condition n\ = N/2 we obtain 



2ze 



0n _ 



N 
M' 



(8) 



since the number of leftover monomers is n\ = zMe^ M . 
Combining this equation with Eq. (J7J yields the critical 
bonding strength 

N r- 

/3Ag cb =^—- In z - ln(2 - y/2). (9) 

For our simulations, N/M = 0.0385, so Ag cb = -4.5 k B T 
if we take z = 6. 

The theory also allows us to estimate the bonding 
strength required for the mass formation of p-segment 
chains. When \Ag\ increases from 0, starts to 
grow, eventually surpasses m, and becomes the domi- 
nant mode. Then n% decreases and gets passed by n^, 
and likewise for n p with larger p. Therefore n p (p > 2) 
shows a maximum at certain Ag c (p) that depends on p. 
We identify this Ag c (p) as the critical bonding strength 
that signals the assembly transition to p-segment chains. 
The corresponding equation is 



drip 
dA~g 



= 0. 



(10) 



Combining this condition with Eqs. ((5$ and (O yields 



p-1 



P 



1' 



(11) 



and a simple expression for the critical bonding energy 



PAg c (p) =ln 



N 
M 



In z — In 



1 



(12) 



Note that |Ap c (2)| = 4.76 k B T is larger than |A3 cb | that 
signals the overall assembly transition. Also note that 
Ag c {p) — > — oo whenp — > 1, which seems counterintuitive 
but is indeed the byproduct of the approximate F in 
Eq. (Ql from the lattice model. 

The above theory is developed for the straight poly- 
merization of wedge monomers through vertical bonding. 
Our simulations further indicate that it also works in the 
case of ring formation, i.e., the case where Al ^ and 
Ay = 0. The main difference is that a ring only contains 
certain number of monomers. In this case, p m ax is finite. 
However, it turns out that in the range of Al where rings 
do form, the results from a finite p ma x, such as p m ax = 13 
for ideal rings, and those from an infinite p max are very 
close. The reason is that chains with p > 13 contribute 
little to the summation in Eq. (JSJ, which is the constraint 
used to determine n p . 

To further describe the formation of tubes, we need to 
take into account simultaneously both the vertical (chain 
formation) and lateral (ring formation) bonding interac- 
tions. It turns out a simple picture can reasonably ex- 
plain our simulation results as presented in the next sec- 
tion. The basic idea is that as long as both the vertical 



and lateral bonding can occur and stay stable (its mean- 
ing will be clarified in the next section), then it is possi- 
ble for monomers to self-assemble into tubes, through a 
cluster growth process. For example, for a system domi- 
nated by p-segment linear chains formed through vertical 
bonding, the tubes will emerge as long as the side-to-side 
(lateral) bonding of these chains can occur and stay sta- 
ble. In this case, the required lateral bonding strength 
will depend on p, in addition to the p-dependence of the 
vertical bonding Ag c (p) that controls the transition to 
p-segment chains. With a mapping between Ag and A, 
this will lead to a (Av,Al) curve that signals the tran- 
sition to the tubular structure phase from the filament 
phase. In another example, if Al is large enough to in- 
duce rings, then tubes will appear when Ay is strong 
enough to lead to stable bonding that stacks rings to- 
gether. As shown later, these arguments lead to accurate 
estimates of critical bonding strengths that are required 
for the tube formation and they agree well with our MD 
calculations. 



IV. RESULTS 
A. Structure Diagram 

We have varied the interaction strengths Al and Ay 
primarily to determine the appropriate values to form 
tubules. Given the experience with capsid assembly^ 
we expect that there will be a narrow range of interac- 
tion strengths that are good for forming tubules without 
defects. We are also interested in the competition be- 
tween the lateral and vertical interactions. In natural 
MTs, it was viewed for a long time that long protofila- 
ments formed and merged to form the tubule^ ' 41 ! 42 In 
other words, the vertical interaction was stronger and 
vertical growth dominated. More recent experiments find 
that the growing end of a MT can vary by lengths cor- 
responding to only a few tubulin dimers or even just one 
dimer, which implies that the vertical and lateral growth 
and interaction strengths are not so distinct 42r— In cells, 
MTs typically grow from a preexisting nucleation center, 
e.g., the ring complex involving 7-tubulinj^ There are 
no such preformed nucleation sites in our simulations. 
However, small clusters formed at an earlier stage of the 
assembly serves such a role as long as they can stay sta- 
ble for a substantial interval. During the assembly, such 
clusters capture free monomers or even other clusters to 
grow into tubes or larger clusters. More details will be 
included in the later subsection on assembly kinetics. 

The wedge monomers do form tubules as shown in 
Fig. [5J This image shows the system at Al = 4.4 k B T 
and Ay — 2.6 k B T, which is one of the parameter sets 
for which tubule structures readily form. This system 
is a large one that has 5000 monomers and many, long 
tubules are seen to have formed. One unexpected result 
is that there are multiple tubules that are helical. For 
example, the pink tubule near the lower front of the im- 
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FIG. 2: Image of system with A L = 4.4 k B T and Ay = 
2.6 ksT showing many tubules formed. The wedge monomers 
are represented by spheres. Monomers in the same cluster 
have the same color. There are only 32 colors; some colors 
are used more than once. 



age clearly shows the helical turn of monomers instead of 
the straight stacking of rings. The helical turn can also 
be seen clearly in the cyan tubule just above and to the 
right of the pink tubule. 

Overall, the MD simulations find an interesting set of 
different structures as shown in Fig. [3] Images of the cor- 
responding structures are shown in Fig. 01 The images 
are colored based on cluster size. Two wedge monomers 
are considered to be in the same cluster if their mutual 
bonding energy is larger than 50% of the full bonding 
energy 4 A. We have analyzed the bonding energy dis- 
tributions and this criterion has accurately captured all 
bonds and does not introduce any artificial ones. 

In the lower-left region of the structure diagram the 
interaction strengths are too small for assembly to oc- 
cur. The majority of the system is in the single monomer 
state. A snapshot of such a system is given in Fig. Ufa) 
(M state). Along the Al = line, standard polymer- 
ization of a straight living polymer occurs, once the 
value of Ay is large enough to yield stable bonds. For 
Ay > 3.2 k B T, stable oligomers start to emerge but sin- 
gle monomers still dominate. We denote this state as 
Oy but consider it only as a subset of the M state. At 
Ay = 4.4 k B T the number of single monomers decreases 
to 50% of the total number of the starting monomers. 
Similarly, along the Ay — line, oligomers appear when 
Al > 3.2 k B T and at Al > 4.4 kBT the monomer count 
is below 50%. Between A L = 4.4 k B T and 5 k B T, the 
polymers are partial rings or arcs (A state, see Fig.[4|b); 
the interaction strength is not large enough for full rings. 
Beyond Al = 5 k B T, individual rings form (R state, see 
Fig. He)). 

Moving above the Ay — line for A L > 3.0 k B T, the 
structure does not change until a large enough value is 
reached to achieve vertical assembly. For Al values that 




FIG. 3: The structure diagram of the assembly of the wedge- 
shaped monomers from MD simulations (symbols) and com- 
parison with the prediction of the thermodynamic theory 
(lines): monomers (M) (open diamonds), oligomers (Ol and 
Ov) (solid circles), partial rings or arcs (A) (downward open 
triangles), full rings (R) (open circles), full tubules (T) (solid 
squares), filaments (F) (upward open triangles), large sheets 
(S) (open squares), multiple clusters (C) (pluses), kinetically 
trapped percolated cluster (K) (crosses). 



have produced even partial rings, the combination of two 
partial rings occurs at low values of Ay , since the multi- 
ple monomers bind between the partial rings yielding the 
extra energy needed to form a larger cluster. Full tubules 
occur in the region denoted with solid squares in Fig. [3] 
(T state) and an example is shown in Fig. 0|d). With 
respect to the ease of assembly, we find that tubules are 
more well formed (i.e. fewer defects) in the region circled 
by the ellipse in the figure, where Al > Ay. The long 
tubules (blue and green) in Fig. H{d) are helical, which 
was not expected given that the wedge was designed to 
form nonhelical tubes. 

At relatively small Al, in the region where Ay > 
4 k B T, oligomers and long filaments form as vertical 
bonding dominates. One example of this filament state 
(labeled F) is shown in Fig. gfe). When A L > 1 k B T, fil- 
aments can bond side-by-side and curved sheets form (S 
state, see Fig. SJf)). O nce the value of Ay > 5 ~ 6 k B T, 
sheets grow very long due to very strong vertical bond- 
ing, and while rings of monomers do form within some 
of the curved sheets as long as Al > 2 k B T, well formed 
tubules do not appear. At even larger values of both 
Al and Ay, the system tends to get kinetically trapped. 
Most or even all monomers join a dominating cluster that 
percolates the whole simulation box (K state, see Fig. 
S£h)); a diffusion limited aggregation is most likely oc- 
curring in this region. The K state is identified by the 
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(e) (f) (g) (h) 



FIG. 4: Images of states in structure diagram: (a) the monomer state (M state) at Al = Av = 3.0 k B T; (b) partial rings and 
arcs (A state) at Al = 4.4 k B T and Av = 1.7 k B T; (c) rings (R state) at Al = 5.5 k B T and Ay — 0.6 k B T; (d) assembly of 
tubes (T state) at A L = 3.9 k B T and Av = 2.6 k B T; (e) filaments (F state) at A L = 1.1 k B T and = 5.2 fc B T; (f) large 
sheets (S state) at A L = 1.6 k B T and = 5.2 k B T. (g) multiple clusters (C state) at A L = 4.4 k B T and Av = 3.3 k B T. 
(h) kinetically trapped percolated cluster state (K state) at Al = 4.4 k B T and Av = 5.2 /cbT. For clarity, each monomer is 
represented as a sphere, and free monomers are only shown except in the panel (a). Monomers are colored based on the size 
of the cluster they reside in (red to blue). Note that periodic boundary conditions are in effect in all three directions so that 
monomers/clusters on opposite sides of the box might actually belong to the same cluster. 



criterion that the largest cluster contains more than 50% 
of the all monomers. When A L is kept large and Ay is 
lowered to around 4 /cbT, systems are still in the kinet- 
ically trapped phase. However, now there exist multiple 
smaller clusters in the system (C state, see Fig. Big)). 
These clusters are very stable and have many defects in- 
duced by the strong lateral bonding. It is very difficult 
for them to fit to each other and fuse into tubes. 

It must be emphasized that there are no sharp bound- 
aries amongst the S, K, C, and T states. In the cross-over 
region between these states, very often various structures 
can coexist with each other. For example, in the C state, 
there are several clusters, inside which closed rings may 
appear and there may also exist tubes. However, these 
tubes usually contain a lot of defects, which can stay for 
a long time. Furthermore, the distinction between the 
C state and the K state is just that in the K state the 
dynamics results in clusters that form early in time, and 
then collide and merge into a percolated structure. 

Large Ay does not favor the formation of tubes be- 
cause long filaments form first, and when they initially 
bind, they are not aligned parallel. While they may in 
time align parallel, the filaments tend to be offset, and the 
process for two filaments to diffuse to remove the offset is 
very slow at best and not observed in the simulation. In 
this case, even though the tube phase has the lowest free 



energy and is the thermodynamically equilibrium state, 
it is not achieved at least on the MD timescale. 



We also observed that when Al is very large and Ay 
is small, there emerges interesting structures other than 
rings. One example is shown in Fig. [SJ which is a long 
helical twist formed at A L = 8.8 k-^T and Ay = 0. 
The number of monomers in this kind of twist can be 
much larger than 13. Furthermore, at a large Al where 
rings can be formed, the number of monomers in a turn 
varies substantially. For example, at Al — 5.5 kBT and 
Ay = 0.1 ksT, we have found rings containing 12 to 14 
monomers with a peak at 13. At Al = 6.6 and 
Ay = 0.1 kBT the number can vary anywhere from f f to 
16. The rings with a large number of monomers usually 
have a slightly twisted shape out of a plane because of the 
geometrical constraint of the wedge size and shape. In 
these cases, tubes do form when Ay is increased. How- 
ever, the tubes tend to have many defects, partly because 
the strong bonding interactions make the relaxation very 
difficult, which hinders the removal of structural defects 
through adjustment of bonds. 
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FIG. 5: A long helical twist formed at Al = 8.8 ksT and 
A v = 0. 



B. Comparison with Thermodynamic Theory 

A basic issue is the minimal well depth (A) necessary to 
achieve assembly. The structure diagram (Fig. [3]) shows 
that assembly does not occur until A^y is greater than 
about 3 k B T . Given that the monomers are rigid bodies 
with 6 total degrees of freedom, the entropy per monomer 
is 3 k B T. Thus, we do not expect dimer formation until 
beyond this range of A. How large a magnitude of A 
must be to form dimers that are stable enough to grow 
into trimers is an open question, which we have addressed 
for our system. To determine the binding energy to form 
dimers we examined in more detail the interaction be- 
tween a pair of wedges. These dimer calculations provide 
a mapping between the model parameters Alt/ and the 
bonding energy Ag in the thermodynamic theory out- 
lined in Sec. IIIII At the level of dimer formation, simula- 
tions show there is no essential difference between lateral 
and vertical bonding. We directly calculated the bonding- 
energy between two monomers starting in a bound dimer 
state, bonded either laterally or vertically, respectively. 
Below we use A to designate either Al or Ay. Simula- 
tion run times were from 5 x 10 3 t to 5 x 10 4 t; the longer 
run times were for cases near or beyond the stability limit 
where the bond had longer lifetimes. 
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FIG. 6: The bonding energy C/(f) as a function of time for 
A — 3 ksT. The horizontal line shows the mean bonding 
energy U b before the dissociation, and the time tb is the bond 
lifetime in this MD run. 



noted as U(t) to indicate its temporal dependence. A 
typical example of U(t) is shown in Fig.|H]for A = 3 k B T. 
The debonding transition to the unbound state is quite 
sharp, which allows us to define a bonding time tb that 
describes the life time of the bond. The mean bonding 
energy from the plateau of U(t) at t < tb is defined as 
the energy per bond U b , which is shown as the horizontal 
line in Fig. O 

By starting with the same initial dimer configuration 
and following different paths in the phase space (i.e. vary- 
ing random number seed in Langevin thermostat), we 
have calculated the mean bonding time (tb) as a function 
of the soft potential strength A. Results on (tb) vs. A 
are shown in the main panel of Fig. [TJa). The inset of 
Fig. EJa) shows ln(tf,) vs. A, from which a change in the 
slope (from 1.05 to 2.61) can be easily identified when A 
increases beyond about 2.6 k B T. A corresponding expo- 
nential fit to (tb) ~ A is also shown in the main panel 
and indicates that the bond life time (tb) starts to grow 
faster with A at A/k B T > 2.6. 

Figure [S] shows that the mean energy U b per bond is 
less than the full bonding value 4A because of thermal 
fluctuations. More data on Ub at various A's are shown 
in Fig. [TJb). In addition to calculating Ub from dimer 
simulations, we also did calculations for a starting state 
of a ring of 13 monomers bonding laterally as a check 
for many-body effects and found none. The energy per 
bond did not depend on the starting state as seen in 
Fig. Efb). Data in Fig. [T](b) indicates a best fit Ub — 
4* A — 3.41 k B T. The slope has the expected value 4 that 
is independent of temperature T. The negative intercept 
—3.41 k B T reflects the effect of thermal motions. 

To establish the mapping between A and Ag, we cal- 
culate the chain distributions n p of p-segment chains for 
systems along the line Al — 0, i.e., the straight polymer- 
ization case. Along this line, dimers and short oligomers 
start to appear when Ay > 3 k B T, but longer chains 
will not form significantly until Ay > 4 ^ 5 k B T. The 
distributions of n p are calculated in MD simulations at 
various Ay 's and compared with predictions of the ther- 
modynamic theory (Eq.fJI)])); in the latter the depen- 
dence on Ag enters in. The results for pn p vs. p at 
Ay = 4.4, 4.8, and 5.2 k B T are shown in Fig. [8] Chain 
distributions shown here are for systems starting with 
all monomers free. To verify that they actually represent 
the equilibrium distribution at the corresponding Ay , we 
did another type of simulations where the starting state 
was many pre-formed tubes. For the values of Ay used 
here, these tubes started to disassemble into filaments, 
oligomers, and monomers with time. It was confirmed 
that these systems reached the same final distribution of 
Tip as those reported here starting with free monomers. 

The theoretical lines from Eq. with various Ag's 
are also shown in Fig. [5] to compare with the MD results. 
We found that the agreement is quite satisfactory if we 
use the the relation 



The instantaneous bonding energy of the dimer is de- 



-Ag = U B - 9.62 k B T 

= 4*A- 13.03 k B T, 



(13) 
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FIG. 7: (a) The mean bonding time (tt) vs. A. Data are cal- 
culated with dimers starting at full bonding either laterally 
(circles) or vertically (squares). The inset is the linear- log plot 
of the same data set in the main panel. The dotted (solid) line 
in the inset is the linear fit to the first (last) 6 data points. The 
solid line in the main panel is the corresponding exponential 
curve for the last 6 data points, (b) The bonding potential 
Ub vs. A for one bond calculated with various starting config- 
urations and temperatures: two monomers bonding laterally 
at T = LOe/fce (circles); two monomers bonding vertically at 
T — 1.0e/fcB (squares); two monomers bonding laterally at 
T — 0.1e/fcB (triangles). Both Ub and A are normalized by 
k B T. The solid line shows a linear fit Ub = 4 * A - 3.41 k B T. 



which effectively introduces a minimal energy of 9.62 k^T 
for the bonding interactions of two monomers to form a 
dimer. With this mapping, the thermodynamic theory 
accurately describes the straight polymerization process. 
It is also found that assembly along the line Ay = 
is well described by the thermodynamic theory. Further 
note that as Ay is increased, the maximum of n p shifts to 
larger p, which validates the criterion (Eq. (jlpp ) we em- 
ployed earlier to determine the critical strength, Ag c (p), 
for the transition to p-segment chains. 

The physical origin of the minimal energy ~ 13 fee? 1 f° r 
stable bonding incorporates three factors. For any bond 
to occur, the bonding energy must be stronger than the 
entropic cost which is 3 k^T. Even for bonding energies 
above this, there is a time in which the dimer will dis- 
sociate (see Fig. [7]). For growth of a chain, let alone a 
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FIG. 8: The distribution of p-segment chains for Al = 
and various Ay's: Av = 4.4 feT (circles); Ay = 4.8 k B T 
(squares); Ay = 5.2 k B T (triangles). Lines are calculated 
from Eq. © and Eq. © with Ag = -(4 * Ay - 13.03 k B T). 



tubule, the stability time must be greater than at least 
the collision time, which is the minimal time required 
for an additional monomer to bind. The collision time 
depends on the concentration. Thus, the value of A de- 
termined here is only valid for the single concentration 
treated here. Similar results have been found in the cal- 
culation for a simpler patchy particle model by Sciortino 
et al.£^ Their work treats a hard core sphere with two 
short-ranged attractive spots at the particle poles, which 
can be treated analytically. Figure 4 in Ref. 36 shows 
the polymerization transition as a function of the well 
depth and density. For this simple model, the calcula- 
tions can be done for a wide range of densities. The value 
of the well depth necessary for polymerization depends 
logarithmically on the density. At a density of 0.01, the 
minimum well depth for polymerization is about 10 k^T. 
In our simulations, along Al = polymerization occurs 
at about A = 3.3 fceT or a (surface-surface) well depth 
of 13.2 /cbT. The two systems are geometrically different 
in multiple ways and should not produce identical val- 
ues. However, both produce values that are many kBT 
beyond 3 kBT, which is related to the time scales of dif- 
fusion and dissociation. The translational diffusion times 
should be very similar for the two systems. Because our 
wedge monomers have an orientation dependence, the ro- 
tational diffusion times for binding will be longer than the 
sticky sphere model, which will result in a larger value 
for the well depth as we find. Our flat surface also causes 
more loss of entropy upon binding than for two spheres, 
which after binding still have rotational degrees of free- 
dom. The binding energy for the wedge system will have 
to have an additional energy to compensate for this rel- 
ative loss of the entropy. Overall, the large well depth 
found in our wedge system is consistent with expecta- 
tions and the magnitude is reasonable for such a system. 
Furthermore, Erickson also estimated the work required 
to immobilize a subunit in a dimer or polymer and ob- 
tained a range from 11.7 to 18.4 k B T^- within which the 
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minimal energy calculated here is located. 

With the relationship between A and Ag, we have 
calculated boundaries between the different structures 
found in Fig. [3] using the thermodynamic theory with 
some simple but ad hoc arguments^ In the thermody- 
namic theory we obtained the critical binding energy 
Ag c (p) that determines the assembly transition to p- 
segment chains. From Ag c (p), we can easily compute the 
corresponding soft-potential strength A c at the transition 
through Ag c (p) = -(4 * A c - 13.03 k B T) (for simplicity 
the energy unit k B T is not always explicitly included in 
the following discussion). 

When Al is small, no lateral bonding occurs and the 
problem reduces to a living polymerization phenomenon. 
As discussed in Sec. IIIIl the onset of assembly is iden- 
tified as when half of the monomers end up in the as- 
sembled state. This leads to a critical bonding en- 
ergy Ag c b = —4.5 k B T, corresponding a soft potential 
strength A cb = 4.4 k B T. For A v < 4.4 k B T, most of the 
system is in the single monomer state. Once Ay > 4.4 
k^T, significant self-assembly occurs with most of the 
monomers now in some n-mer. Similarly, when Ay is 
small there is no vertical bonding and significant self- 
assembly of monomers into partial rings through lateral 
bonding occurs when A^ > 4.4 k B T. In Fig. [3] the lines 
at A L = 4.4 k B T and at A v = 4.4 k B T are from this cal- 
culation and they correspond well with the boundaries 
between the oligomer state and A state, and between the 
oligomer state and the F state as determined from the 
MD simulations. 

In the region corresponding to the single monomer M 
state (including the oligomer states Ol and Oy), there 
is a sub-region bounded by the dashed-line parts of the 
two curves that separate the possible tube state from 
other states. Inside this sub-region, monomers can only 
form unstable dimers from time to time and the number 
of such dimers is small at any instant. The boundaries 
of this region roughly corresponds to Ag < 0, which 
from Eq. (fTS")) is A < 3.3 k B T. These boundaries are 
also consistent with the earlier estimate of A for a stable 
bond from the calculation of bond life time t\, as shown in 
Fig. E£a), which shows t& becomes nonzero near 3 k B T. 
A better estimate will be calculated below from the as- 
sembly kinetics and it also indicates that 3.3 k B T is the 
boundary value. 

Full ring formation determines the boundary between 
the A and R states. Since an ideal ring constrains p to be 
13, we can use the condition n p (p = 13) = 1 to estimate 
the lower bound of this boundary. This gives us a critical 
soft potential strength Al — 5 k B T which is shown as a 
line in Fig. [3J and fits the simulations data for boundary 
with the R state. 

At Ag — Ag c (p), the p-segment chains (or partial 
rings) dominate. Then in order for these chains to 
assemble into tubes, the bonding potential that holds 
two chains side by side has to be larger than 9.62/p 
since there are p side-to-side bonds, and this leads to 
A > (3.41 + 9.62/p)/4, where everything is in the unit of 



k B T. The curves given by Ag c (p) = -(4*A 1 - 13.03) 
and A 2 > (3.41 + 9.62/p)/4 with (A 1 ,A 2 ) = (A L , Ay) or 
(Ai,A 2 ) = (Av,Al) determines the boundary between 
the F and S states and between the A/R states and the T 
state, respectively. These are the curved lines in Fig. [3] 

As shown in Fig. [3J there is a remarkable agreement 
between the structure diagram calculated from MD sim- 
ulations and the prediction of the thermodynamic the- 
ory. The predicated boundaries at M, F, S, A, R, and T 
states all fit the simulation results well. This indicates 
that the main feature of the assembly behavior of wedge- 
shaped monomers is captured by the simple free energy 
developed in Sec. IIIIl Figure [3] further shows that the 
kinetically trapped states (the S, C, K states) dominate 
when Ay is large. The thermodynamic theory can not 
capture this noncquilibrium behavior. 



C. Assembly Kinetics 

We have shown that wedge-shaped monomers can 
self-assemble into various structures, depending on the 
strengths of the lateral and vertical bonding interaction. 
When we start with a system with all monomers un- 
bonded, the number of free monomers decreases over time 
with the progress of self-assembly. The rate of consump- 
tion of free monomers is controlled by the strength of 
the bonding potential. In Fig. [^a) , results for the num- 
ber of free monomers n\ normalized by the total num- 
ber of monomers (500 here) is plotted against time t for 
A L = Ay = A, where A ranges from 3.3 to 4.4 k B T. The 
strength 3.3 k B T is close to the threshold for assembly, 
and the corresponding n\ decreases very slowly on this 
time scale, though at a later time it will eventually drops 
below 50% of starting free monomers. At A = 3.4 k B T, 
n\ clearly decreases with time but many free monomers 
remain at the end of the time range shown here. As 
A increases, the rate of reduction becomes larger. At 
A = 4.4 k B T, ni decreases to over a period of only 
3 x 10 5 r, which is an order of magnitude shorter than 
the typical time of MD runs used to calculate the struc- 
ture diagram. 

Beyond an initial drop in n\ at very short times, the n\ 
vs. t curves can be fit to a functional form a*exp(— t/t{), 
with a and t\ as the fitting parameters. The time decay 
constant for the single monomer count is t\. This time 
constant is for data beyond a very short initial time, when 
n\{t) drops sharply. From Fig. [H^a) it is clear that even 
for the A — i.ik B T case the limit to t = is n\ < 1. 
The fits are shown as smooth lines in Fig. [DJa). The plot 
of the time constant t\ vs. A in Fig. [9^b) shows that t\ 
appears to reach a critical value near A ~ 3.2 k B T, which 
is consistent with the earlier estimate of 3.3 k B T as the 
threshold for some dimers to form. 
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FIG. 9: (Color Online) (a) The fraction of free monomers vs. 
time for Al = Av = A. For the first eight curves from top to 
bottom, A increases from 3.3 k&T to 4.0 feT in increments of 
0.1 fceT. The bottommost curve is for A — 4.4e. The smooth 
lines are the fits to a * exp(— t/ii). (b) The fitting parameter 
ti vs. A. The line is a guide to the eye. 



D. Helicity of Assembled Tubes 

The visualization of the formation of an assembled 
structure reveals important aspect of the assembly pro- 
cess. The consecutive snapshots of assembly of one 
tubule in the system for Fig. [5] shown in Fig. [TU1 illustrate 
an interesting example. In Fig. UOf a). two partial rings 
stacking vertically have already formed, and a dimer is 
approaching this cluster. In Fig.fTUTb). the dimer has col- 
lided and bound to the partial double ring. In Fig.fTUtc). 
the binding of the dimer results in the formation of a full 
ring. Then in Fig. HUT d). one wedge in the closed ring 
oscillates noticeably and eventually in Fig. ITUf e). this 
wedge separates from its neighbor and slips to bond with 
the wedges in the partial ring below. In this manner 
a helical ring is formed. These helical structures typi- 
cally have 12 or 13 monomers per turn. This particular 
structure grows into a larger helical tubule with several 
complete turns before the simulation ends. 

Another interesting assembly dynamics involves the 
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FIG. 10: Consecutive snapshots at Al = 4.4 &bT and 
Ay = 2.6 ksT show the formation of a closed ring ((a)-(c)) 
by capturing a dimer and the subsequent relaxation ((d)-(e)) 
into a helical structure. 
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FIG. 11: Consecutive snapshots show the capture of a cluster 
by a partial tube and the formation of a closed one ((a)- 
(f)). The tube subsequently relaxes into a helical structure 
((g)-(i)). In (j), each monomer is represented by a sphere to 
illustrate clearly the helical feature of the tube. 



formation of tubes not by the addition of monomer and 
dimers, but by collisions between large clusters. Fig- 
ure [Tl] shows an example of this dynamics; the two clus- 
ters are displayed through two different shades of gray. In 
Fig. [TTT a). the larger cluster contains two stacked rings 
(nonhclical) in the top half and 3 stacked arcs on the 
bottom half. This cluster is colliding with a curved sheet 
in this first image. In Fig. ITTT b) the curved sheet ro- 
tates and orients to fit into the missing bottom part of 
the larger cluster. In Fig. [TTTc) the two clusters begin to 
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merge. Note that the right hand side of the curved sheet 
(dark cluster) is straight as this point. In Fig. [TTTd') . 
the two clusters have merged and undergone further re- 
arrangement. In particular, the wedge monomers on the 
right hand side, which in Fig. ITTT c) are straight, on top 
of each other, have now translated in a staircase fash- 
ion. Figure lllf e) shows that the triangular gap above 
one of the dark wedges has healed and the whole com- 
bined structure is now a helical tube. The final image in 
Fig. mT f) shows the wedges as spheres, clearly displaying 
a well formed helical tube. 

The geometry of our wedges are designed to form a 
13-monomer ring and then rings can stack to form tubes. 
Though we do observe nonhelical tubes with 13 filaments, 
we also see nonhelical tubes with only 12 monomers in 
each of their rings. However, helical tubes containing 12 
or 13 filaments are more frequently found in our simu- 
lations. We have directly compared the total potential 
energy of two preformed tubes of the same number of 
monomers with one helical and the other not. The non- 
helical tube does have the lower energy. However, the 
difference in energy is small and about equal to the fluc- 
tuation in energy. Thus, helical structure is sufficiently 
close in energy that it is not surprising that we see so 
many. 



V. DISCUSSION 

With dynamics that involves assembly and disassem- 
bly of structures, achieving equilibrium in a simulation 
is typically a difficult task. This is true in some parts of 
the structure diagram, but the boundaries between the 
structures have been well determined at least to the pre- 
cision that we desire. The decay in the number of single 
monomers (Fig. [5]) shows that the self-assembly process 
is very slow near the onset of assembly. We can not run 
the simulations at A = 3.3 k-gT for long enough times 
to achieve equilibrium distributions. However, by exam- 
ining the kinetics as in Fig. [5Jb) , we know that there 
is a boundary near 3.3 k^T. The data presented here 
provides sufficient information for anyone interested in 
more details of a particular structure to know what pa- 
rameter range to focus on, which is our interest. We 
have achieved equilibrium distributions for some parts of 
the structure diagram. In particular, the distributions in 
Fig. [5] are equilibrium values. We demonstrated equi- 
libration by performing simulations of assembly (start- 
ing from monomers) and of disassembly (starting from 
tubules) . 

At large values of A, kinetics tends to be the driver. 
For the kinetically trapped region, the main point is that 
this is a region to avoid when trying to form tubules. 
Even experimentally, where much longer time scales can 
be easily reached, there is no advantage to work in this 
region, because the kinetically trapped structures become 
trapped at relatively short times. 

Experimentally, the advantage of longer equilibration 



times is used coincidentally with lower concentrations. 
In this work, we have kept the concentration constant in 
order to keep the number of variables manageable. The 
choice of the particular concentration was a compromise 
between lower concentrations which would yield better 
assembly and higher concentrations at which assembly 
progresses faster. It is apparent that the structure dia- 
gram will shift when concentration is varied. However, 
the thermodynamic theory predicts that the binding en- 
ergy Ag only depends on concentration logarithmically. 36 
Thus the qualitative picture of various structures pre- 
sented here should still remain valid, with some minor 
quantitative changes in the bond strengths required for 
tubule assembly. 

Our results indicate that tubes are more easily and bet- 
ter formed with few defects when the lateral binding is 
slightly stronger than the vertical binding, i.e., Al > Ay. 
The desired values of Al is around 4 k^T and Ay around 
3 fceT, corresponding to bonding energies about 13 k^T 
laterally and 9 ksT vertically, respectively. These num- 
bers should serve as reasonable guides for experimental 
design of wedge-shaped molecules that can self-assemble 
into tubular structures. We have also found that very 
long tubes can form when Ay = 3.9 fceT and Al = 2.6 
k^T . In this case, only a few tubes and many monomers 
coexist almost without intermediate clusters. Tubes con- 
taining over 300 monomers and having greater than 20 
turns have been assembled in the simulations. However, 
tubes grow very slowly in this case. 

There has been significant interest in modeling MT 
assembly kinetics and dynamics 4 ^ - — , including treating 
mechanical aspect a 54 ' 55 . In Ref. 51, VanBuren et al. es- 
timated bond energies within the MT lattice. They pre- 
dicted that the standard free energy gain of the forma- 
tion of a vertical bond is 6.8 to 9.4 fceT including the 
entropy contribution of immobilizing a dimer in the MT 
lattice. Using our estimate of entropy loss of a bond 
formation as 9.6 k^T, this corresponds to an energy of 
16.4 to 19 ksT for a vertical bond, which is translated to 
Ay — 5.0 ~ 5.6 fceT in our model. They also estimated 
that a lateral bond is much weaker with an energy of 3.2 
to 5.7 ksT, i.e., Al — 1.7 ~ 2.3 k^T. In our structure 
diagram (Fig. [3]), the region encircled by these values 
of Al and Ay accommodates large sheets. Though the 
thermodynamic theory predicts that the lowest energy 
state is for these sheets to close into tubes, such clo- 
sure is hindered by the quick formation of many sheets 
and depletion of free monomers, and difficulty of two or 
more sheets to join together to fuse into tubes because of 
structural mismatch. Only on a time scale much longer 
than the one achievable with MD, can bond breaking and 
structural adjustment occur to allow tubes to form. Also 
note that VanBuren et al. assumed that all the free en- 
ergy of immobilizing a dimer is provided by the vertical 
bonding. If instead we assumed that this free energy cost 
is jointly compensated by both the vertical and the lat- 
eral bonds, then the lateral bond energy (i.e., Al) would 
become larger while the vertical one (i.e., Ay) smaller, 
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which moves the estimate of bond energies of a MT lat- 
tice closer to the region where we observe tube formation 
in our simulations. 

In a recent stochastic simulation work of Wu et al. 
that studied the intermediate sheet structure found dur- 
ing MT assembly, two types of lateral bonds and one 
type of vertical bonds were assumed to have energies of 
13 to 17.5 k B T and 19 k B T, respectively^ These cor- 
respond to A L = 4.1 ~ 5.2 k B T and Ay = 5.6 k B T. 
The lateral bond strength is inside the range found in 
this paper where tubular structures are formed, but the 
vertical bond strength is in the region where kinetically 
trapped structures are found in our simulations. How- 
ever, the monomer concentration in their work is an or- 
der of magnitude lower than ours, which implies that a 
larger value of Ay is required for tube formation in their 
simulations because of the increase of diffusion time and 
the decrease of collision rate. Furthermore, in stochas- 
tic simulations it is much easier for a system to escape 
from kinetic traps and reach a thermodynamic equilib- 
rium state where tubes are expected to form. We thus 
conclude that the parameter set used by Wu et al. is 
consistent with the range of bond strengths suitable for 
tube formation found in this paper. 

Besides the interaction strengths, the geometry of 
the monomer is also important in determining the self- 
assembled structure. Surprisingly, we have found that 
helical tubes are more common than nonhelical tubes. 
The nonhelical tube has the lowest energy, but the dif- 
ference between the two states is small, which is found 
less than the fluctuation in the energy. Thus helical struc- 
tures emerge with entropy as a driven force. The forma- 
tion of the helical tubes occurs in a variety of ways as 
shown in Figs. ITUl and [TO Once the helical structure has 
formed for a length of about one turn, then it is stable 
as helical and grows from the ends. The collisional refor- 
mation seen in Fig. [TT]is not a common mechanism. The 
most common mechanism appears to be the formation of 
a curved sheet of a few layers that has sufficient filaments 
to bend and touch itself. In most cases the self-contact 
is with a twist yielding helical tubes. This mechanisms 
has been seen in other systems such as lipid nanotubes. 

We mentioned that when Al > 5 k B T and Ay < 1 k B T 
rings with various number (11 ~ 16) of monomers can 
form. In our model, 13 monomers fit to make a per- 
fect ring. So rings containing more than 13 monomers 
have to take a somewhat out-of-plane shape because our 
monomer is a rigid composite body. This shape makes 
tube formation difficult, though tubes do form when Ay 
is increased above 1 k B T in this range of Al. How- 
ever, only tubes with 12 or 13 filaments, either helical or 
nonhelical, are observed; the majority has 13 filaments. 
This is slightly different from the case of MT assembly 
in vitro, where the protofilament number of MTs varies 
between 10 and 15, with the vast majority having 14 
protofilaments^ The smaller range found in our simula- 
tions is clearly due to the fact that our monomer is rigid. 
The distortion of the rings containing other than 12 or 13 



monomers is incompatible with the geometry of a tube, 
and the monomer cannot adjust its shape to reduce the 
distortion of the ring. 

We have used the wedge shape to promote formation 
of the tubule structure. We believe that this is not neces- 
sary. A more spherical or ellipsoidal shape could be used 
with the lateral interaction sites being located off center 
in order to promote ring formation. However, changing 
to a more curved shape would have some consequences 
which would require adjustments. On the one hand, 
the flat wedge surfaces yield a tighter fit between bound 
monomers, which raises the entropic cost and must be 
compensated for in the binding energy. On the other 
hand, a more flexible bond will be more unstable (lower 
barrier) and potentially allow other composite geome- 
tries. We have seen in the helical twist that even small 
variations in the geometry between two bound monomers 
can yield very different structures on larger scales. The 
wedge shape brings some simplification which enables 
such difficult simulations, but some of the future work 
will need to match more closely the actual molecular 
building blocks. 

The goal of this work is to study a model monomer 
designed to form tubules, which is motivated by an in- 
terest in developing artificial MTs. This model can be 
further developed to study real MTs and it is interesting 
to see how well the model compares to MT assembly and 
structure, though only a few aspects of tubulin dimers 
are incorporated in the model, which will limit a direct 
comparison. The basic result of forming tubules and find- 
ing the appropriate range of interaction strengths, which 
was discussed above, is the necessary first step. The 
fact that helical tubules are formed in the simulations 
is a pleasant surprise, although one that is understand- 
able. Furthermore, it indicates that some features in the 
monomer must be adjusted in order to control the helical 
assembly. The preference for the lateral strength to be 
greater than the vertical strength is interesting and de- 
serves more study. The results show the importance of 
nucleation in MT assembly, but generally speaking this 
is not a surprise in biological systems which tend to be 
controlled in a systematic manner. 



VI. CONCLUSIONS 

We have developed a minimal model of tube formation 
from wedge-shaped monomers. Our simulations have 
identified a rich structure diagram of the assembly, which 
agrees well with the prediction of a simple Flory-Huggins 
lattice theory. Results indicate that tubes form when 
both the attractive lateral and vertical binding interac- 
tion energies are at the scale of ~ 10 k B T. Tubes with 
fewer defects are more easily formed when the lateral in- 
teraction is slightly stronger than the vertical one. The 
stacking of rings is less susceptible to kinetic traps than 
the alignment and piecing together of filaments. There is 
more than one assembly mechanism for tube formation. 
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We have observed that in addition to the main growth 
mechanism by absorbing free monomers, tubes can also 
grow through capturing dimers and oligomers, or even 
through collisions between clusters. Our results also re- 
veal that helical tubes are more frequently formed than 
nonhelical tubes, despite the fact that our monomer ge- 
ometry and interaction potentials favor nonhelical struc- 
tures, which have the lowest energy. Part of the reason 
for this is that the energy difference is small, which allows 
transitions between rings and helices of a single turn, for 
example. 

This work sets the stage for further model develop- 
ments that will incorporate additional features necessary 
for simulating other aspects of MTs such as strong de- 
polymerization. Our interest is more in developing design 
rules for constructing artificial MTs, but that involves un- 
derstanding MTs in a more general manner. The present 
results indicate the best range of interaction strengths 
to achieve artificial tubule formation. In addition, the 
simulations show that cells have to control multiple as- 
pects of the formation dynamics to prevent defects from 
forming in the assembly of MTs. Tubules can collide and 
bind, which tends to form structures with defects. In the 
helical tubules, we find some with a pitch of 2 monomers 
instead of 1. All these factors suggest that instability 
inherent within MTs is also a mechanism to handle un- 
wanted, poorly growing structures. 



the wedge monomers, z = 6 for a cubic lattice is most 
appropriate. Thus the number of different ways to place 
the first segment of the i-th p-segment chain is 



K x =z [M - M p - (i - l)p] 



(16) 



Since each monomer can form two bonds on its oppo- 
site faces in linear chains, the second segment has to 
be placed in 2 neighboring cells of the first segment 
where bonding is possible. Using mean-field approxi- 
mation, the probability for one such cell to be empty 
is [M — M p — (i — l)p — 1\/M, The number of ways to 
place the second segment of the z-th p-segment chain is 
then 



Ko = 



2 [A/ - M p - (i - l)p - 1] 
M ' 



(17) 



Now the third segment can only be placed in the cell 
adjacent to the one that the second segment resides, 
but opposite to the cell that the first segment re- 
sides. The probability for this cell to be available is 
[M - M p - (i- l)p - 2] /M. The number of ways to 
place the third or any segment j > 3 of the i-th p-segment 
chain is 



M-Mp-{i-l)p -{j-l) 
M 



(18) 



From Eq. (1141) we obtain 



Appendix 

In this appendix, we give a detailed derivation of the 
free energy used in the main text. Since chains of dif- 
ferent lengths are distinguishable, we can first place the 
1-segment chains (leftover monomers) onto the lattice, 
then place the 2-segment chains (dimers), and so on. Let 
£pi be the number of ways to place the i-th p-segment 
chain, and W p be the total number of ways to place all 
p-segment chains, and taking into account the indistin- 
guishability amongst the p-segment chains, we obtain 



p »=i j=i 



i 



z 



(M - M p )\ 



\ M n p (p-l) (M-M p+ i)!' 



(19) 



Here the identity M p+ i = M p + pn p is used. The factor 
1/2 is included because there are 2 ways to choose the 
first segment of a rigid linear chain, and thus there is a 
2-fold degeneracy in the above counting of ways to place 
chains. Finally, the total number of different ways to 
place all chains is 



(14) 



The number of occupied sites by all chains up to but not 
including the p-segments is 



p-i 



M, 



j > 



for p > 2 



(15) 



Note that Mi = 0. There are M — M p cells available 
for the first p-segment chain. When we start placing 
the i-th p-segment chain, there are M — M p — (i — l)p 
cells empty for its first segment (monomer). Since the 
wedge monomer is an anisotropic object, its orientation 
contributes to the entropy of the system. In the lattice 
model, the monomer has z possible orientations, where z 
is the coordination number of the lattice. For modeling 



z=\[w p 

P =i 

Ml 



(20) 



where n s — M — X)p=i x P n p = M — N is the number of 
solvent (empty) cells after all chains are placed. 

The configurational entropy S c of the assembled sys- 
tem is 



S c = fee In Z 

Pmax 

P =i 



(rip 



i ] n p i A ( 21 ) 

in z — n„ m — - + n Tj ' 

p M 



p 



where is the Boltzmann constant, and terms that only 
depend on n s and M are discarded since they are con- 
stants. 
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The Helmholtz free energy can be written as United States Department of Energy under Contract No. 

DE-AC04-94AL85000. This research was supported by 

F = E -TS, 
_ 1 



the U.S. Department of Energy, Office of Basic Energy 
2 Pmax ,<22) Sciences, Division of Materials Sciences and Engineering 
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P =i 



where (p — l)Ag is the total energy gain to form a p- 
segment chain, which has p—1 bonds, and j3 = 1/knT. 
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